function [ar_coeffs,sigma2] = ar_cov(x,p)
N = length(x);
% 协方差法估计 AR 系数
R = zeros(p+1, p+1); % 初始化自相关协方差矩阵
for j = 0:p
    for k = 0:p
        for n = max(j,k)+1:N
            R(j+1,k+1) = R(j+1,k+1) + x(n)*conj(x(n - max(j,k) + min(j,k)));
        end
    end
end

R = R/(N - p); % 归一化
r = R(2:end, 1); % 自相关向量
R = R(2:end, 2:end); % 去除第一行第一列

% R = zeros(p, p);
% r = zeros(p, 1);

% for k = 1:p
%     for l = 1:p
%         for n = max(k, l):N
%             R(k, l) = R(k, l) + conj(x(n - k)) * x(n - l);
%         end
%     end
%     for n = k:N
%         r(k) = r(k) + conj(x(n - k)) * x(n);
%     end
% end

% 求解 Yule - Walker 方程得到 AR 系数
ar_coeffs = -R \ r; 

% 计算预测误差功率
sigma2 = R(1, 1) + ar_coeffs' * r;
% ar_coeffs = transpose(ar_coeffs);
ar_coeffs = [1,transpose(ar_coeffs)]; % same with aryule()

end

